covid project preliminary analysis

knitr::opts_chunk$set(echo = TRUE)

Introduction

In this file, we cleaned and explored our following datasets:
- CDC Social Vulnerability Index/SVI (2018): https://www.atsdr.cdc.gov/placeandhealth/svi/data_documentation_download.html
- JHU covid data (curren): https://github.com/CSSEGISandData/COVID-19
- Google mobility reports (current): https://www.google.com/covid19/mobility/?hl=en
- NCSL unemployment data: https://www.ncsl.org/research/labor-and-employment/state-unemployment-update.aspx
The CDC data provided us with demographic information such as poverty rates, incomes, and minority rates. Even though the SVI surveys were conducted The JHU covid date provided us with historical covid progressions. Mobility and unemployment data both provided us relevant information during the pandamic. All datasets are analyzed on state-level and had been geographically categorized into Midwest, Northeast, South, and West based on convention/Wikipedia.

For convenice, jump to Graphs to view visual representation.

Data cleaning

Importing packages

library(tidyverse) 
library(dplyr)     

Data cleaning for CDC data

selecting relevant indices (refer to codebook)

svi.raw <- read_csv("./datasets/SVI2018_US_COUNTY.csv")
svi.raw %>% 
  select(STATE,
         COUNTY,
         FIPS,
         AREA_SQMI,
         E_TOTPOP,
         E_HU,
         E_HH,
         E_POV,
         E_UNEMP,
         E_PCI,
         E_NOHSDP,
         E_AGE65,
         E_DISABL,
         E_MINRTY,
         E_LIMENG,
         E_MOBILE,
         E_CROWD,
         E_NOVEH,
         E_GROUPQ,
         E_UNINSUR) -> svi.data

-999 means unavailable data so replaced with 0

svi.selected <- svi.data
svi.selected[svi.selected == "-999"] <- NA
# only New Mexico has NA responses for est. PCI, poverty, and unemployment

further cleaning svi data and calculating total amount of variables

svi.selected %>% 
  # PCI is per capita income 
  # calculating total income to derive statewide per capita income later 
  mutate(income = E_TOTPOP * E_PCI) %>% 
  group_by(STATE) %>% 
  summarise(area = sum(AREA_SQMI),
            # aggregating to state level
            total.population = sum(E_TOTPOP),
            total.housing.units = sum(E_HU),
            total.household = sum(E_HH),
            total.poverty = sum(E_POV),
            total.unemployed = sum(E_UNEMP),
            total.income = sum(income),
            # dividing into statewide per capita
            total.capita = total.income/total.population,
            total.no.HS = sum(E_NOHSDP),
            total.age65 = sum(E_AGE65),
            total.disabled = sum(E_DISABL),
            total.minority = sum(E_MINRTY),
            total.eng.proficiency = sum(E_LIMENG),
            total.mobile = sum(E_MOBILE),
            total.crowd = sum(E_CROWD),
            total.no.vehicle = sum(E_NOVEH),
            total.institutionalized = sum(E_GROUPQ),
            total.uninsured = sum(E_UNINSUR)) -> svi.state

further cleaning svi data and calculating percentages of variables

svi.state %>%
  mutate(population.density = total.population/area,
         poverty.rate = total.poverty/total.population,
         unemployed.rate = total.unemployed/total.population,
         per.capita = total.capita,
         no.HS.rate = total.no.HS/total.population,
         age65.rate = total.age65/total.population,
         disabled.rate = total.disabled/total.population,
         minority.rate = total.minority/total.population,
         poor.English.rate = total.eng.proficiency/total.population,
         mobile.home.rate = total.mobile/total.population,
         household.crowd.rate = total.crowd/total.population,
         no.vehicle.rate = total.no.vehicle/total.population,
         institutionalized.rate = total.institutionalized/total.population,
         uninsured.rate = total.uninsured/total.population) -> svi.state

Data cleaning for unemployment data

unemployment.2020.raw <- read_csv("./datasets/unemployment2020.csv")
unemployment.2020.data <- 
  unemployment.2020.raw[,colSums(is.na(unemployment.2020.raw))<nrow(unemployment.2020.raw)]
unemployment.2020 <- na.omit(unemployment.2020.data)
svi.state %>%
  select(STATE, unemployed.rate) -> unemployment.2018

svi.state %>%
  select(-unemployed.rate) -> svi.state

svi.state$STATE <- str_to_title(svi.state$STATE[1:51])

# write.csv(svi.state,"C:\\Users\\fan\\Desktop\\preliminary cdc data version 1.csv")

combining unemployment data

unemployment <- cbind(unemployment.2018, unemployment.2020) %>% 
  mutate(unemployed.rate = round(unemployed.rate * 100, digit = 1)) %>% 
  rename(unemployment.rate.2018 = unemployed.rate) %>% 
  select(-c(STATE)) %>% # select everything but 
  rename(state = State,
         Jan = Jan.,
         Feb = Feb.,
         Mar = March,
         Apr = April,
         May = May,
         Jun = June,
         Jul = July,
         Aug = Aug.,
         Sep = Sept.) %>%
  pivot_longer(-c(state, unemployment.rate.2018),
               names_to = "month",
               values_to = "unemployment.rate.2020") %>% 
  mutate(month = match(month, month.abb))

unemployments <- unemployment[c(2, 3, 1, 4)]

# write.csv(unemployments,"E:/UCSD/2020-2021 Senior/PSYC 201 Project/CDC + Employment Data/unemployment rate version 1.csv")

JHU covid data cleaning

loading JHU covid death historical dateset

covid.death.raw <- read_csv("./datasets/csse_covid_19_time_series/time_series_covid19_deaths_US.csv")
# glimpse(covid.death.raw)

cleaning JHU covid death historical dateset

# Exclude these states
# American Samoa, Diamond Princess, Grand Princess, Guam, Northern Mariana Islands, Puerto Rico, Virgin Islands
covid.death <- covid.death.raw %>% 
  select(-c(UID, iso2, iso3, code3, FIPS, Admin2, Country_Region, Lat, Long_, Combined_Key)) %>%
  rename(state = Province_State, population = Population) %>% 
  group_by(state) %>% # Specify group indicator
  summarise(across(everything(), sum)) %>% # Specify column
  pivot_longer(-c(state, population),
               names_to = "date",
               values_to = "cumulative.death.cases") %>% 
  filter(state %in% c(state.name, "District of Columbia")) %>% 
  group_by(state) %>% 
  mutate(death.cases = cumulative.death.cases - lag(cumulative.death.cases),
         date = as.Date(date, format = "%m/%d/%y"),
         month = match(months(date), month.name),
         death.cases = replace_na(death.cases, 0))

# fix negative cumulative cases
while (any(covid.death$death.cases < 0)) {
  covid.death <- covid.death %>% 
    mutate(cumulative.death.cases = ifelse(cumulative.death.cases < lag(cumulative.death.cases), 
                                           lag(cumulative.death.cases), # if smaller than previous
                                           cumulative.death.cases), # if not
           cumulative.death.cases = replace_na(cumulative.death.cases, 0),
           death.cases = cumulative.death.cases - lag(cumulative.death.cases),
           death.cases = replace_na(death.cases, 0))
}

covid.deaths <- covid.death[c(1, 2, 3, 6, 4, 5)]

loading JHU covid confirmed historical dateset

covid.confirmed.raw <- read_csv("./datasets/csse_covid_19_time_series/time_series_covid19_confirmed_US.csv")
# glimpse(covid.confirmed.raw)

cleaning JHU covid confirmed historical dateset

covid.confirmed <- covid.confirmed.raw %>% 
  select(-c(UID, iso2, iso3, code3, FIPS, Admin2, Country_Region, Lat, Long_, Combined_Key)) %>%
  rename(state = Province_State) %>% 
  group_by(state) %>% # Specify group indicator
  summarise(across(everything(), sum)) %>% # Specify column
  pivot_longer(-c(state),
               names_to = "date",
               values_to = "cumulative.confirmed.cases") %>% 
  filter(state %in% c(state.name, "District of Columbia")) %>% 
  group_by(state) %>% 
  mutate(confirmed.cases = cumulative.confirmed.cases - lag(cumulative.confirmed.cases),
         date = as.Date(date, format = "%m/%d/%y"), 
         confirmed.cases = replace_na(confirmed.cases, 0))

# fix negative cumulative cases
while (any(covid.confirmed$confirmed.cases < 0)) {
  covid.confirmed <- covid.confirmed %>% 
    mutate(cumulative.confirmed.cases = ifelse(cumulative.confirmed.cases < lag(cumulative.confirmed.cases), 
                                         lag(cumulative.confirmed.cases), # if smaller than previous
                                         cumulative.confirmed.cases), # if not
           cumulative.confirmed.cases = replace_na(cumulative.confirmed.cases, 0),
           confirmed.cases = cumulative.confirmed.cases - lag(cumulative.confirmed.cases),
           confirmed.cases = replace_na(confirmed.cases, 0))
}

# combining death and confirmed into one dataset
covid <- covid.deaths %>% inner_join(covid.confirmed, by = c("state","date"))

Looping through dates

dates <- format(seq(from = as.Date("2020/04/13"), to = as.Date("2020/10/24"), by = "day"), "%m-%d-%Y")
date = "04-12-2020"
reports.raw <- read_csv(paste("./datasets/csse_covid_19_daily_reports_us/", date, ".csv", sep = ""))
report <- reports.raw %>%
  select(-c(Country_Region, Last_Update, Lat, Long_, Confirmed, Deaths, FIPS, UID, ISO3)) %>% 
  rename(state = Province_State,
         recovered.cases = Recovered,
         active.cases = Active, 
         incident.rate = Incident_Rate, # positive results per 100,000 persons.
         people.tested = People_Tested,
         people.hospitalized = People_Hospitalized,
         mortality.rate = Mortality_Rate,
         testing.rate = Testing_Rate,
         hospitalization.rate = Hospitalization_Rate) %>% 
  filter(state %in% c(state.name, "District of Columbia")) %>% 
  mutate(date = as.Date(date, format = "%m-%d-%y"),
         # month = match(months(date), month.name),
         incident.rate = incident.rate/1000,
         # based on documentation (100k population)
         testing.rate = testing.rate/1000) %>%
    select(-c(people.hospitalized, hospitalization.rate, recovered.cases))
  


for (date in dates){
  reports.raw <- read_csv(paste("./datasets/csse_covid_19_daily_reports_us/", date, ".csv", sep = ""))
  rpt <- reports.raw %>%
    select(-c(Country_Region, Last_Update, Lat, Long_, Confirmed, Deaths, FIPS, UID, ISO3)) %>% 
    rename(state = Province_State,
           recovered.cases = Recovered,
           active.cases = Active,
           incident.rate = Incident_Rate,
           people.tested = People_Tested,
           people.hospitalized = People_Hospitalized,
           mortality.rate = Mortality_Rate,
           testing.rate = Testing_Rate,
           hospitalization.rate = Hospitalization_Rate) %>% 
    filter(state %in% c(state.name, "District of Columbia")) %>% 
    mutate(date = as.Date(date, format = "%m-%d-%y"),
           # month = match(months(date), month.name),
           incident.rate = incident.rate/1000,
           testing.rate = testing.rate/1000) %>%
    select(-c(people.hospitalized, hospitalization.rate, recovered.cases))
  ### North Dekota has > 100% testing rate 
  report <- rbind(report, rpt)
}

reports <- report[c(1, 6, 7, 2, 3, 4, 5)][order(report$state, report$date),]

covid.time.series <- covid %>% 
  full_join(reports, by = c("state","date")) %>% 
  mutate(active.rates = 100*(active.cases / population)) # added normalized active rates

# covid.time.series <- reports
# View(covid.time.series)
#write.csv(covid.time.series, "E:/UCSD/2020-2021 Senior/PSYC 201 Project/CJHU Databse/covid time series version 1.csv")

Cleaning mobility data

Importing Google mobility dataset

US_Region_2020_mobility <- read.csv( "./datasets/2020_US_Region_Mobility_Report.csv", header=T, na.strings=c("","NA"))

Cleaning mobility dataset

US_State_2020_mobility<- US_Region_2020_mobility %>% filter(sub_region_1 != "Null")
US_State_2020_mobility_without_county <- US_State_2020_mobility %>% filter(iso_3166_2_code != "Null")

# unique(US_State_2020_mobility_without_county$sub_region_1)
US_State_2020_mobility_without_county$date <- as.Date(US_State_2020_mobility_without_county$date, "%m/%d/%y")


month <- as.numeric(format(US_State_2020_mobility_without_county$date, "%m"))

StateMonthMobility <- US_State_2020_mobility_without_county %>% mutate(month = month)


StateMonthmobility <- StateMonthMobility %>% filter(retail_and_recreation_percent_change_from_baseline != "NA")%>% 
  filter(grocery_and_pharmacy_percent_change_from_baseline != "NA")%>%
  filter(parks_percent_change_from_baseline != "NA")%>%
  filter(transit_stations_percent_change_from_baseline!= "NA")%>% 
  filter(workplaces_percent_change_from_baseline!= "NA")%>% 
  filter(residential_percent_change_from_baseline!= "NA")

StateMonthMobility2 <- StateMonthMobility %>% 
  group_by(month,sub_region_1) %>% 
  summarise(mean_retial_and_recreation_percent_change = 
              mean(retail_and_recreation_percent_change_from_baseline),
            mean_grocery_and_pharmacy_percent_change = 
              mean(grocery_and_pharmacy_percent_change_from_baseline), 
            mean_parks_percent_change = mean(parks_percent_change_from_baseline),
            mean_transit_stations_percent_change = 
              mean(transit_stations_percent_change_from_baseline),
            mean_workplaces_percent_change = 
              mean(workplaces_percent_change_from_baseline),
            mean_residential_percent_change= 
              mean(residential_percent_change_from_baseline))

Visualization codes

Geographic categorization

Midwest.states <- c("North Dakota","South Dakota","Nebraska","Kansas","Indiana","Missouri","Iowa","Ohio","Wisconsin","Michigan","Minnesota","Illinois")
Northeast.states <- c("Pennsylvania","New Hampshire","Maine","Connecticut","New Jersey","Rhode Island","New York","Vermont","Massachusetts")
South.states <- c("West Virginia","Oklahoma","Kentucky","Alabama","Arkansas","Tennessee","Louisiana","Mississippi","South Carolina","Texas","Georgia","North Carolina","Florida","Virginia","Delaware","Maryland","District of Columbia")
West.states <- c("Wyoming","Idaho","Montana","Utah","Alaska","Arizona","Nevada","Colorado","New Mexico","Oregon","Washington","California","Hawaii")

JHU covid data

Geographic categories for covid data based on Wikipedia

# Categorizing states to plot in 4 sub graphs
covid.state.midwest <- filter(covid.time.series, state %in% Midwest.states)
covid.state.northeast <- filter(covid.time.series, state %in% Northeast.states)
covid.state.south <- filter(covid.time.series, state %in% South.states)
covid.state.west <- filter(covid.time.series, state %in% West.states)

covid death graphs

covid.death.midwest <- covid.state.midwest %>%
  ggplot(aes(x = date, y = cumulative.death.cases, color = state)) + 
  geom_line() +
  scale_color_hue(name='Midwest States')+
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(0,40000,5000), limits=c(0,40000),
                     labels=as.character(seq(0,40000,5000)))+
  ylab("Cumulative Death Cases") +
  xlab("Date") + 
  ggtitle("Cumulative Death Cases in Midwest US") + 
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style
  
covid.death.northeast <- covid.state.northeast %>%
  ggplot(aes(x = date, y = cumulative.death.cases, color = state)) + 
  geom_line() +
  scale_color_hue(name='Northeast States')+
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(0,40000,5000), limits=c(0,40000),
                     labels=as.character(seq(0,40000,5000)))+
  ylab("Cumulative Death Cases") +
  xlab("Date") + 
  ggtitle("Cumulative Death Cases in Northeast US") + 
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

covid.death.south<- covid.state.south %>%
  ggplot(aes(x = date, y = cumulative.death.cases, color = state)) + 
  geom_line() +
  scale_color_hue(name='South States')+
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(0,40000,5000), limits=c(0,40000),
                     labels=as.character(seq(0,40000,5000)))+
  ylab("Cumulative Death Cases") +
  xlab("Date") + 
  ggtitle("Cumulative Death Cases in South US") + 
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

covid.death.west <- covid.state.west %>%
  ggplot(aes(x = date, y = cumulative.death.cases, color = state)) + 
  geom_line() +
  scale_color_hue(name='West States')+
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(0,40000,5000), limits=c(0,40000),
                     labels=as.character(seq(0,40000,5000)))+
  ylab("Cumulative Death Cases") +
  xlab("Date") + 
  ggtitle("Cumulative Death Cases in West US") + 
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

covid confirmed cases graphs

covid.confirm.midwest <- covid.state.midwest %>%
  ggplot(aes(x = date, y = cumulative.confirmed.cases, color = state)) + 
  geom_line() +
  scale_color_hue(name='Midwest States')+
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  ylim(0, 1000000)+
  ylab("Cumulative Confirmed Cases") +
  xlab("Date") + 
  ggtitle("Cumulative Confirmed Cases in Midwest US") + 
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

covid.confirm.northeast <- covid.state.northeast %>%
  ggplot(aes(x = date, y = cumulative.confirmed.cases, color = state)) + 
  geom_line() +
  scale_color_hue(name='Northeast States')+
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  ylim(0, 1000000)+
  ylab("Cumulative Confirmed Cases") +
  xlab("Date") + 
  ggtitle("Cumulative Confirmed Cases in Northeast US") + 
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

covid.confirm.south <- covid.state.south %>%
  ggplot(aes(x = date, y = cumulative.confirmed.cases, color = state)) + 
  geom_line() +
  scale_color_hue(name='South States')+
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  ylim(0, 1000000)+
  ylab("Cumulative Confirmed Cases") +
  xlab("Date") + 
  ggtitle("Cumulative Confirmed Cases in South US") + 
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

covid.confirm.west <- covid.state.west %>%
  ggplot(aes(x = date, y = cumulative.confirmed.cases, color = state)) + 
  geom_line() +
  scale_color_hue(name='West States')+
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  ylim(0, 1000000)+
  ylab("Cumulative Confirmed Cases") +
  xlab("Date") + 
  ggtitle("Cumulative Confirmed Cases in West US") + 
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

covid testing rate graphs

# Only starting from April 12th because we only have data dtarting from then
covid.test.rate.midwest <- covid.state.midwest %>%
  filter(date > as.Date("2020-04-11")) %>% 
  ggplot(aes(x = date, y = testing.rate, color = state)) + 
  geom_line() +
  scale_color_hue(name='Midwest States')+
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(0,120,20), limits=c(0,105),
                     labels=as.character(seq(0,120,20)))+
  ylab("Testing Rate") +
  xlab("Date") + 
  ggtitle("Testing Rate in Midwest US") + 
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

covid.test.rate.northeast <- covid.state.northeast %>%
  filter(date > as.Date("2020-04-11")) %>% 
  ggplot(aes(x = date, y = testing.rate, color = state)) + 
  geom_line() +
  scale_color_hue(name='Northeast States')+
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(0,120,20), limits=c(0,105),
                     labels=as.character(seq(0,120,20)))+
  ylab("Testing Rate") +
  xlab("Date") + 
  ggtitle("Testing Rate in Northeast US") + 
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

covid.test.rate.south <- covid.state.south %>%
  filter(date > as.Date("2020-04-11")) %>% 
  ggplot(aes(x = date, y = testing.rate, color = state)) + 
  geom_line() +
  scale_color_hue(name='South States')+
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(0,120,20), limits=c(0,105),
                     labels=as.character(seq(0,120,20)))+
  ylab("Testing Rate") +
  xlab("Date") + 
  ggtitle("Testing Rate in South US") + 
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

covid.test.rate.west <- covid.state.west %>%
  filter(date > as.Date("2020-04-11")) %>% 
  ggplot(aes(x = date, y = testing.rate, color = state)) + 
  geom_line() +
  scale_color_hue(name='West States')+
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(0,120,20), limits=c(0,105),
                     labels=as.character(seq(0,120,20)))+
  ylab("Testing Rate") +
  xlab("Date") + 
  ggtitle("Testing Rate in West US") + 
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

covid testing postive rate graphs

# Only starting from April 12th because we only have data starting from then
covid.incident.rate.midwest <- covid.state.midwest %>%
  filter(date > as.Date("2020-04-11")) %>% 
  ggplot(aes(x = date, y = incident.rate, color = state)) + 
  geom_line() +
  scale_color_hue(name='Midwest States')+
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  ylim(0, 5)+
  ylab("Incident Rate") +
  xlab("Date") + 
  ggtitle("Incident Rate in Midwest US") + 
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

covid.incident.rate.northeast <- covid.state.northeast %>%
  filter(date > as.Date("2020-04-11")) %>% 
  ggplot(aes(x = date, y = incident.rate, color = state)) + 
  geom_line() +
  scale_color_hue(name='Northeast States')+
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  ylim(0, 5)+
  ylab("Incident Rate") +
  xlab("Date") + 
  ggtitle("Incident Rate in Northeast US") + 
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

covid.incident.rate.south <- covid.state.south %>%
  filter(date > as.Date("2020-04-11")) %>% 
  ggplot(aes(x = date, y = incident.rate, color = state)) + 
  geom_line() +
  scale_color_hue(name='South States')+
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  ylim(0, 5)+
  ylab("Incident Rate") +
  xlab("Date") + 
  ggtitle("Incident Rate in South US") + 
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

covid.incident.rate.west <- covid.state.west %>%
  filter(date > as.Date("2020-04-11")) %>% 
  ggplot(aes(x = date, y = incident.rate, color = state)) + 
  geom_line() +
  scale_color_hue(name='West States')+
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  ylim(0, 5)+
  ylab("Incident Rate") +
  xlab("Date") + 
  ggtitle("Incident Rate in West US") + 
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

Google mobility data

Geographic categories for covid data based on Wikipedia

StateMonthMobility_Midwest <- filter(StateMonthMobility2,sub_region_1 %in% Midwest.states)
StateMonthMobility_Northeast <- filter(StateMonthMobility2,sub_region_1 %in% Northeast.states)
StateMonthMobility_South <- filter(StateMonthMobility2,sub_region_1 %in% South.states)
StateMonthMobility_West <- filter(StateMonthMobility2,sub_region_1 %in% West.states)

Plotting for Mean Percent Change of Retail+Recreation

mobility.r.midwest <- StateMonthMobility_Midwest %>% 
  mutate(month = as.Date(paste("2020-", month,"-01", sep=""))) %>% 
  ggplot(aes(y=mean_retial_and_recreation_percent_change, 
             x=month, color=sub_region_1)) +
  geom_point() + 
  geom_line() + 
  geom_path()+
  ggtitle("Mean Percent Change of Retail and Recreation in Midwest US")+
  ylab("Mean Percent Change of Retail and Recreation") + 
  xlab("Month") + 
  labs(color="Midwest States") + 
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(-60,60,10), limits=c(-65,20),
                     labels=as.character(seq(-60,60,10)))+
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

mobility.r.northeast <- StateMonthMobility_Northeast %>% 
  mutate(month = as.Date(paste("2020-", month,"-01", sep=""))) %>% 
  ggplot(aes(y=mean_retial_and_recreation_percent_change, 
             x=month, color=sub_region_1)) +
  geom_point() + 
  geom_line() + 
  geom_path()+
  ggtitle("Mean Percent Change of Retail and Recreation in Northeast US")+
  ylab("Mean Percent Change of Retail and Recreation") + 
  xlab("Month") + 
  labs(color="Northeast States") + 
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(-60,60,10), limits=c(-65,20),
                     labels=as.character(seq(-60,60,10)))+
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

mobility.r.south <- StateMonthMobility_South %>% 
  mutate(month = as.Date(paste("2020-", month,"-01", sep=""))) %>% 
  ggplot(aes(y=mean_retial_and_recreation_percent_change, 
             x=month, color=sub_region_1)) +
  geom_point() + 
  geom_line() + 
  geom_path()+
  ggtitle("Mean Percent Change of Retail and Recreation in South US")+
  ylab("Mean Percent Change of Retail and Recreation") + 
  xlab("Month") + 
  labs(color="South States") + 
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(-60,60,10), limits=c(-65,20),
                     labels=as.character(seq(-60,60,10)))+
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style


mobility.r.west <- StateMonthMobility_West %>% 
  mutate(month = as.Date(paste("2020-", month,"-01", sep=""))) %>% 
  ggplot(aes(y=mean_retial_and_recreation_percent_change, 
             x=month, color=sub_region_1)) +
  geom_point() + 
  geom_line() + 
  geom_path()+
  ggtitle("Mean Percent Change of Retail and Recreation in West US")+
  ylab("Mean Percent Change of Retail and Recreation") + 
  xlab("Month") + 
  labs(color="West States") + 
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(-60,60,10), limits=c(-65,20),
                     labels=as.character(seq(-60,60,10)))+
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

# grid.arrange(p1,p2,p3,p4,nrow=2,top="Relationship between Mean Percent Change of Retail+Recreation and Month by States")

Plotting for Mean Percent Change of Grocery+Pharmacy

mobility.gp.midwest <- StateMonthMobility_Midwest %>% 
  mutate(month = as.Date(paste("2020-", month,"-01", sep=""))) %>% 
  ggplot(aes(y=mean_grocery_and_pharmacy_percent_change, 
             x=month, color=sub_region_1)) +
  geom_point() + 
  geom_line() + 
  geom_path()+
  ggtitle("Mean Percent Change of Grocery and Pharmacy in Midwest US")+
  ylab("Mean Percent Change of Grocery and Pharmacy") + 
  xlab("Month") + 
  labs(color="Midwest States") + 
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(-60,60,10), limits=c(-35,35),
                     labels=as.character(seq(-60,60,10)))+
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

mobility.gp.northeast <- StateMonthMobility_Northeast %>% 
  mutate(month = as.Date(paste("2020-", month,"-01", sep=""))) %>% 
  ggplot(aes(y=mean_grocery_and_pharmacy_percent_change, 
             x=month, color=sub_region_1)) +
  geom_point() + 
  geom_line() + 
  geom_path()+
  ggtitle("Mean Percent Change of Grocery and Pharmacy in Northeast US")+
  ylab("Mean Percent Change of Grocery and Pharmacy") + 
  xlab("Month") + 
  labs(color="Northeast States") + 
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(-60,60,10), limits=c(-35,35),
                     labels=as.character(seq(-60,60,10)))+
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

mobility.gp.south <- StateMonthMobility_South %>% 
  mutate(month = as.Date(paste("2020-", month,"-01", sep=""))) %>% 
  ggplot(aes(y=mean_grocery_and_pharmacy_percent_change, 
             x=month, color=sub_region_1)) +
  geom_point() + 
  geom_line() + 
  geom_path()+
  ggtitle("Mean Percent Change of Grocery and Pharmacy in South US")+
  ylab("Mean Percent Change of Grocery and Pharmacy") + 
  xlab("Month") + 
  labs(color="South States") + 
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(-60,60,10), limits=c(-35,35),
                     labels=as.character(seq(-60,60,10)))+
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

mobility.gp.west <- StateMonthMobility_West %>% 
  mutate(month = as.Date(paste("2020-", month,"-01", sep=""))) %>% 
  ggplot(aes(y=mean_grocery_and_pharmacy_percent_change, 
             x=month, color=sub_region_1)) +
  geom_point() + 
  geom_line() + 
  geom_path()+
  ggtitle("Mean Percent Change of Grocery and Pharmacy in West US")+
  ylab("Mean Percent Change of Grocery and Pharmacy") + 
  xlab("Month") + 
  labs(color="West States") + 
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(-60,60,10), limits=c(-35,35),
                     labels=as.character(seq(-60,60,10)))+
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

# grid.arrange(q1,q2,q3,q4,nrow=2,top="Relationship between Mean Percent Change of Grocery+Pharmacy and Month by States")

Plotting for Mean Percent Change of Parks

mobility.parks.midwest <- StateMonthMobility_Midwest %>% 
  mutate(month = as.Date(paste("2020-", month,"-01", sep=""))) %>% 
  ggplot(aes(y=mean_parks_percent_change, 
             x=month, color=sub_region_1)) +
  geom_point() + 
  geom_line() + 
  geom_path()+
    ggtitle("Mean Percent Change of Parks in Midwest US")+
  ylab("Mean Percent Change of Parks") + 
  xlab("Month") + 
  labs(color="Midwest States") + 
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(-100,400,50), limits=c(-75,400),
                     labels=as.character(seq(-100,400,50)))+
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

mobility.parks.northeast <- StateMonthMobility_Northeast %>% 
  mutate(month = as.Date(paste("2020-", month,"-01", sep=""))) %>% 
  ggplot(aes(y=mean_parks_percent_change, 
             x=month, color=sub_region_1)) +
  geom_point() + 
  geom_line() + 
  geom_path()+
    ggtitle("Mean Percent Change of Parks in Northeast US")+
  ylab("Mean Percent Change of Parks") + 
  xlab("Month") + 
  labs(color="Northeast States") + 
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(-100,400,50), limits=c(-75,400),
                     labels=as.character(seq(-100,400,50)))+
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

mobility.parks.south <- StateMonthMobility_South %>% 
  mutate(month = as.Date(paste("2020-", month,"-01", sep=""))) %>% 
  ggplot(aes(y=mean_parks_percent_change, 
             x=month, color=sub_region_1)) +
  geom_point() + 
  geom_line() + 
  geom_path()+
    ggtitle("Mean Percent Change of Parks in South US")+
  ylab("Mean Percent Change of Parks") + 
  xlab("Month") + 
  labs(color="South States") + 
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(-100,400,50), limits=c(-75,400),
                     labels=as.character(seq(-100,400,50)))+
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

mobility.parks.west <- StateMonthMobility_West %>% 
  mutate(month = as.Date(paste("2020-", month,"-01", sep=""))) %>% 
  ggplot(aes(y=mean_parks_percent_change, 
             x=month, color=sub_region_1)) +
  geom_point() + 
  geom_line() + 
  geom_path()+
    ggtitle("Mean Percent Change of Parks in West US")+
  ylab("Mean Percent Change of Parks") + 
  xlab("Month") + 
  labs(color="West States") + 
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(-100,400,50), limits=c(-75,400),
                     labels=as.character(seq(-100,400,50)))+
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

# grid.arrange(r1,r2,r3,r4,nrow=2,top="Relationship between Mean Percent Change of Parks and Month by States")

Plotting for Mean Percent Change of Transit Stations

mobility.transit.midwest <- StateMonthMobility_Midwest %>% 
  mutate(month = as.Date(paste("2020-", month,"-01", sep=""))) %>% 
  ggplot(aes(y=mean_transit_stations_percent_change, 
             x=month, color=sub_region_1)) +
  geom_point() + 
  geom_line() + 
  geom_path()+
  ggtitle("Mean Percent Change of Transit Stations in Midwest US")+
  ylab("Mean Percent Change of Transit Stations") + 
  xlab("Month") + 
  labs(color="Midwest States") + 
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(-100,100,10), limits=c(-70,45),
                     labels=as.character(seq(-100,100,10)))+
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

mobility.transit.northeast <- StateMonthMobility_Northeast %>% 
  mutate(month = as.Date(paste("2020-", month,"-01", sep=""))) %>% 
  ggplot(aes(y=mean_transit_stations_percent_change, 
             x=month, color=sub_region_1)) +
  geom_point() + 
  geom_line() + 
  geom_path()+
  ggtitle("Mean Percent Change of Transit Stations in Northeast US")+
  ylab("Mean Percent Change of Transit Stations") + 
  xlab("Month") + 
  labs(color="Northeast States") + 
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(-100,100,10), limits=c(-70,45),
                     labels=as.character(seq(-100,100,10)))+
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

mobility.transit.south <- StateMonthMobility_South %>% 
  mutate(month = as.Date(paste("2020-", month,"-01", sep=""))) %>% 
  ggplot(aes(y=mean_transit_stations_percent_change, 
             x=month, color=sub_region_1)) +
  geom_point() + 
  geom_line() + 
  geom_path()+
  ggtitle("Mean Percent Change of Transit Stations in South US")+
  ylab("Mean Percent Change of Transit Stations") + 
  xlab("Month") + 
  labs(color="South States") + 
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(-100,100,10), limits=c(-70,45),
                     labels=as.character(seq(-100,100,10)))+
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

mobility.transit.west <- StateMonthMobility_West %>% 
  mutate(month = as.Date(paste("2020-", month,"-01", sep=""))) %>% 
  ggplot(aes(y=mean_transit_stations_percent_change, 
             x=month, color=sub_region_1)) +
  geom_point() + 
  geom_line() + 
  geom_path()+
  ggtitle("Mean Percent Change of Transit Stations in West US")+
  ylab("Mean Percent Change of Transit Stations") + 
  xlab("Month") + 
  labs(color="West States") + 
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(-100,100,10), limits=c(-70,45),
                     labels=as.character(seq(-100,100,10)))+
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

# grid.arrange(s1,s2,s3,s4,nrow=2,top="Relationship between Mean Percent Change of Transit Stations and Month by States")

Plotting for Mean Percent Change of Workplaces

mobility.work.midwest <- StateMonthMobility_Midwest %>% 
  mutate(month = as.Date(paste("2020-", month,"-01", sep=""))) %>% 
  ggplot(aes(y=mean_workplaces_percent_change, 
             x=month, color=sub_region_1)) +
  geom_point() + 
  geom_line() + 
  geom_path()+
  ggtitle("Mean Percent Change of Workplaces in Midwest US")+
  ylab("Mean Percent Change of Workplaces") + 
  xlab("Month") + 
  labs(color="Midwest States") + 
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(-60,60,10), limits=c(-65,10),
                     labels=as.character(seq(-60,60,10)))+
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

mobility.work.northeast <- StateMonthMobility_Northeast %>% 
  mutate(month = as.Date(paste("2020-", month,"-01", sep=""))) %>% 
  ggplot(aes(y=mean_workplaces_percent_change, 
             x=month, color=sub_region_1)) +
  geom_point() + 
  geom_line() + 
  geom_path()+
  ggtitle("Mean Percent Change of Workplaces in Northeast US")+
  ylab("Mean Percent Change of Workplaces") + 
  xlab("Month") + 
  labs(color="Northeast States") + 
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(-60,60,10), limits=c(-65,10),
                     labels=as.character(seq(-60,60,10)))+
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

mobility.work.south <- StateMonthMobility_South %>% 
  mutate(month = as.Date(paste("2020-", month,"-01", sep=""))) %>% 
  ggplot(aes(y=mean_workplaces_percent_change, 
             x=month, color=sub_region_1)) +
  geom_point() + 
  geom_line() + 
  geom_path()+
  ggtitle("Mean Percent Change of Workplaces in South US")+
  ylab("Mean Percent Change of Workplaces") + 
  xlab("Month") + 
  labs(color="South States") + 
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(-60,60,10), limits=c(-65,10),
                     labels=as.character(seq(-60,60,10)))+
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

mobility.work.west <- StateMonthMobility_West %>% 
  mutate(month = as.Date(paste("2020-", month,"-01", sep=""))) %>% 
  ggplot(aes(y=mean_workplaces_percent_change, 
             x=month, color=sub_region_1)) +
  geom_point() + 
  geom_line() + 
  geom_path()+
  ggtitle("Mean Percent Change of Workplaces in West US")+
  ylab("Mean Percent Change of Workplaces") + 
  xlab("Month") + 
  labs(color="West States") + 
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(-60,60,10), limits=c(-65,10),
                     labels=as.character(seq(-60,60,10)))+
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

# grid.arrange(t1,t2,t3,t4,nrow=2,top="Relationship between Mean Percent Change of Workplaces and Month by States")

Plotting Mean Percent Change of Residential (%change of people staying at home)

mobility.res.midwest <- StateMonthMobility_Midwest %>% 
  mutate(month = as.Date(paste("2020-", month,"-01", sep=""))) %>% 
  ggplot(aes(y=mean_residential_percent_change, 
             x=month, color=sub_region_1)) +
  geom_point() + 
  geom_line() + 
  geom_path()+
  ggtitle("Mean Percent Change of Residential in Midwest US")+
  ylab("Mean Percent Change of Residential") + 
  xlab("Month") + 
  labs(color="Midwest States") + 
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(-60,60,5), limits=c(-5,20),
                     labels=as.character(seq(-60,60,5)))+
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

mobility.res.northeast <- StateMonthMobility_Northeast %>% 
  mutate(month = as.Date(paste("2020-", month,"-01", sep=""))) %>% 
  ggplot(aes(y=mean_residential_percent_change, 
             x=month, color=sub_region_1)) +
  geom_point() + 
  geom_line() + 
  geom_path()+
  ggtitle("Mean Percent Change of Residential in Northeast US")+
  ylab("Mean Percent Change of Residential") + 
  xlab("Month") + 
  labs(color="Northeast States") + 
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(-60,60,5), limits=c(-5,20),
                     labels=as.character(seq(-60,60,5)))+
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

mobility.res.south <- StateMonthMobility_South %>% 
  mutate(month = as.Date(paste("2020-", month,"-01", sep=""))) %>% 
  ggplot(aes(y=mean_residential_percent_change, 
             x=month, color=sub_region_1)) +
  geom_point() + 
  geom_line() + 
  geom_path()+
  ggtitle("Mean Percent Change of Residential in South US")+
  ylab("Mean Percent Change of Residential") + 
  xlab("Month") + 
  labs(color="South States") + 
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(-60,60,5), limits=c(-5,20),
                     labels=as.character(seq(-60,60,5)))+
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

mobility.res.west <- StateMonthMobility_West %>% 
  mutate(month = as.Date(paste("2020-", month,"-01", sep=""))) %>% 
  ggplot(aes(y=mean_residential_percent_change, 
             x=month, color=sub_region_1)) +
  geom_point() + 
  geom_line() + 
  geom_path()+
  ggtitle("Mean Percent Change of Residential in West US")+
  ylab("Mean Percent Change of Residential") + 
  xlab("Month") + 
  labs(color="West States") + 
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(-60,60,5), limits=c(-5,20),
                     labels=as.character(seq(-60,60,5)))+
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

# grid.arrange(u1,u2,u3,u4,nrow=2,top="Relationship between Mean Percent Change of Residential and Month by States")

NCSL unemployment data

Geographic categories for unemployment data based on Wikipedia

unemployment_Midwest <- filter(unemployment, state %in% Midwest.states)
unemployment_Northeast <- filter(unemployment, state %in% Northeast.states)
unemployment_South <- filter(unemployment, state %in% South.states)
unemployment_West <- filter(unemployment, state %in% West.states)

Plotting 2020 unemployment

unemploy.mw <- unemployment_Midwest %>%
  mutate(month = as.Date(paste("2020-", month,"-01", sep=""))) %>% 
  ggplot(aes(x = month, y = unemployment.rate.2020, color = state)) +
  geom_line() +
  geom_point() +
  scale_color_hue(name='Midwest States')+
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(-60,60,5), limits=c(0,30),
                     labels=as.character(seq(-60,60,5)))+
  ylab("Unemployment Rate") +
  xlab("Month") + 
  ggtitle("2020 Unemployment Rate in Midwest US") + 
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

unemploy.ne <- unemployment_Northeast %>%
  mutate(month = as.Date(paste("2020-", month,"-01", sep=""))) %>% 
  ggplot(aes(x = month, y = unemployment.rate.2020, color = state)) + 
  geom_line() +
  geom_point() +
  scale_color_hue(name='Northeast States')+
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(-60,60,5), limits=c(0,30),
                     labels=as.character(seq(-60,60,5)))+
  ylab("Unemployment Rate") +
  xlab("Month") + 
  ggtitle("2020 Unemployment Rate in Northeast US") + 
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

unemploy.s <- unemployment_South %>%
  mutate(month = as.Date(paste("2020-", month,"-01", sep=""))) %>% 
  ggplot(aes(x = month, y = unemployment.rate.2020, color = state)) + 
  geom_line() +
  geom_point() +
  scale_color_hue(name='South States')+
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(-60,60,5), limits=c(0,30),
                     labels=as.character(seq(-60,60,5)))+
  ylab("Unemployment Rate") +
  xlab("Month") + 
  ggtitle("2020 Unemployment Rate in South US") + 
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

unemploy.w <- unemployment_West %>%
  mutate(month = as.Date(paste("2020-", month,"-01", sep=""))) %>% 
  ggplot(aes(x = month, y = unemployment.rate.2020, color = state)) + 
  geom_line() +
  geom_point() +
  scale_color_hue(name='West States')+
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  scale_y_continuous(expand=c(0,1), breaks=seq(-60,60,5), limits=c(0,30),
                     labels=as.character(seq(-60,60,5)))+
  ylab("Unemployment Rate") +
  xlab("Month") + 
  ggtitle("2020 Unemployment Rate in West US") + 
  theme_minimal() +
  theme(axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

CDC SVI data

Geographic categories for cdc data based on Wikipedia

South.states.1 <- c("West Virginia","Oklahoma","Kentucky","Alabama","Arkansas","Tennessee","Louisiana","Mississippi","South Carolina","Texas","Georgia","North Carolina","Florida","Virginia","Delaware","Maryland","District Of Columbia")
# made adjustment for DC spelling difference
svi.state_Midwest <- filter(svi.state, STATE %in% Midwest.states)
svi.state_Northeast <- filter(svi.state, STATE %in% Northeast.states)
svi.state_South <- filter(svi.state, STATE %in% South.states.1)
svi.state_West <- filter(svi.state, STATE %in% West.states)

Plotting poverty rate

svi.pov.mw <- svi.state_Midwest %>%
  ggplot(aes(x = reorder(STATE, poverty.rate), y = poverty.rate, color = STATE, fill = STATE)) + 
  geom_bar(stat = "identity") + 
  coord_flip()+
  ylim(0, 0.21) +
  xlab("State") +
  ylab("Poverty Rate") + 
  ggtitle("Poverty Rate in Midwest US") + 
  theme_minimal() +
  theme(legend.position = 'none', # remove legend,
        axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

svi.pov.ne <- svi.state_Northeast %>%
  ggplot(aes(x = reorder(STATE, poverty.rate), y = poverty.rate, color = STATE, fill = STATE)) + 
  geom_bar(stat = "identity") + 
  coord_flip()+
  ylim(0, 0.21) +
  xlab("State") +
  ylab("Poverty Rate") + 
  ggtitle("Poverty Rate in Northeast US") + 
  theme_minimal() +
  theme(legend.position = 'none', # remove legend,
        axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

svi.pov.s <- svi.state_South %>%
  ggplot(aes(x = reorder(STATE, poverty.rate), y = poverty.rate, color = STATE, fill = STATE)) + 
  geom_bar(stat = "identity") + 
  coord_flip()+
  ylim(0, 0.21) +
  xlab("State") +
  ylab("Poverty Rate") + 
  ggtitle("Poverty Rate in South US") + 
  theme_minimal() +
  theme(legend.position = 'none', # remove legend,
        axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

svi.pov.w <- svi.state_West %>%
  ggplot(aes(x = reorder(STATE, poverty.rate), y = poverty.rate, color = STATE, fill = STATE)) + 
  geom_bar(stat = "identity") + 
  coord_flip()+
  geom_text(x=13, y =0.042, label = "New Mexico data missing",
            size = 5, color = "black")+
  ylim(0, 0.21) +
  xlab("State") +
  ylab("Poverty Rate") + 
  ggtitle("Poverty Rate in South US") + 
  theme_minimal() +
  theme(legend.position = 'none', # remove legend,
        axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style
### New Mexico data missing

Plotting uninsured rate

svi.uninsured.mw <- svi.state_Midwest %>%
  ggplot(aes(x = reorder(STATE, uninsured.rate), y = uninsured.rate, color = STATE, fill = STATE)) + 
  geom_bar(stat = "identity") + 
  coord_flip()+
  ylim(0, 0.175) +
  xlab("State") +
  ylab("Uninsured  Rate") + 
  ggtitle("Uninsured Rate in Midwest US") + 
  theme_minimal() +
  theme(legend.position = 'none', # remove legend,
        axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

svi.uninsured.ne <- svi.state_Northeast %>%
  ggplot(aes(x = reorder(STATE, uninsured.rate), y = uninsured.rate, color = STATE, fill = STATE)) + 
  geom_bar(stat = "identity") + 
  coord_flip()+
  ylim(0, 0.175) +
  xlab("State") +
  ylab("Uninsured  Rate") + 
  ggtitle("Uninsured Rate in Northeast US") + 
  theme_minimal() +
  theme(legend.position = 'none', # remove legend,
        axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

svi.uninsured.s <- svi.state_South %>%
  ggplot(aes(x = reorder(STATE, uninsured.rate), y = uninsured.rate, color = STATE, fill = STATE)) + 
  geom_bar(stat = "identity") + 
  coord_flip()+
  ylim(0, 0.175) +
  xlab("State") +
  ylab("Uninsured  Rate") + 
  ggtitle("Uninsured Rate in South US") + 
  theme_minimal() +
  theme(legend.position = 'none', # remove legend,
        axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

svi.uninsured.w <- svi.state_West %>%
  ggplot(aes(x = reorder(STATE, uninsured.rate), y = uninsured.rate, color = STATE, fill = STATE)) + 
  geom_bar(stat = "identity") + 
  coord_flip()+
  ylim(0, 0.175) +
  xlab("State") +
  ylab("Uninsured  Rate") + 
  ggtitle("Uninsured Rate in West US") + 
  theme_minimal() +
  theme(legend.position = 'none', # remove legend,
        axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

Plotting mobile home rate

svi.mobile.mw <- svi.state_Midwest %>%
  ggplot(aes(x = reorder(STATE, mobile.home.rate), y = mobile.home.rate, color = STATE, fill = STATE)) + 
  geom_bar(stat = "identity") + 
  coord_flip()+
  ylim(0, 0.08) +
  xlab("State") +
  ylab("Mobile Home Rate") + 
  ggtitle("Mobile Home Rate in Midwest US") + 
  theme_minimal() +
  theme(legend.position = 'none', # remove legend,
        axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

svi.mobile.ne <- svi.state_Northeast %>%
  ggplot(aes(x = reorder(STATE, mobile.home.rate), y = mobile.home.rate, color = STATE, fill = STATE)) + 
  geom_bar(stat = "identity") + 
  coord_flip()+
  ylim(0, 0.08) +
  xlab("State") +
  ylab("Mobile Home Rate") + 
  ggtitle("Mobile Home Rate in Northeast US") + 
  theme_minimal() +
  theme(legend.position = 'none', # remove legend,
        axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

svi.mobile.s <- svi.state_South %>%
  ggplot(aes(x = reorder(STATE, mobile.home.rate), y = mobile.home.rate, color = STATE, fill = STATE)) + 
  geom_bar(stat = "identity") + 
  coord_flip()+
  ylim(0, 0.08) +
  xlab("State") +
  ylab("Mobile Home Rate") + 
  ggtitle("Mobile Home Rate in South US") + 
  theme_minimal() +
  theme(legend.position = 'none', # remove legend,
        axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

svi.mobile.w <- svi.state_West %>%
  ggplot(aes(x = reorder(STATE, mobile.home.rate), y = mobile.home.rate, color = STATE, fill = STATE)) + 
  geom_bar(stat = "identity") + 
  coord_flip()+
  ylim(0, 0.08) +
  xlab("State") +
  ylab("Mobile Home Rate") + 
  ggtitle("Mobile Home Rate in West US") + 
  theme_minimal() +
  theme(legend.position = 'none', # remove legend,
        axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

Plotting minority rate

svi.minority.mw <- svi.state_Midwest %>%
  ggplot(aes(x = reorder(STATE, minority.rate), y = minority.rate, color = STATE, fill = STATE)) + 
  geom_bar(stat = "identity") + 
  coord_flip()+
  ylim(0, 0.8) +
  xlab("State") +
  ylab("Minority Rate") + 
  ggtitle("Minority Rate in Midwest US") + 
  theme_minimal() +
  theme(legend.position = 'none', # remove legend,
        axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

svi.minority.ne <- svi.state_Northeast %>%
  ggplot(aes(x = reorder(STATE, minority.rate), y = minority.rate, color = STATE, fill = STATE)) + 
  geom_bar(stat = "identity") + 
  coord_flip()+
  ylim(0, 0.8) +
  xlab("State") +
  ylab("Minority Rate") + 
  ggtitle("Minority Rate in Northeast US") + 
  theme_minimal() +
  theme(legend.position = 'none', # remove legend,
        axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

svi.minority.s <- svi.state_South %>%
  ggplot(aes(x = reorder(STATE, minority.rate), y = minority.rate, color = STATE, fill = STATE)) + 
  geom_bar(stat = "identity") + 
  coord_flip()+
  ylim(0, 0.8) +
  xlab("State") +
  ylab("Minority Rate") + 
  ggtitle("Minority Rate in South US") + 
  theme_minimal() +
  theme(legend.position = 'none', # remove legend,
        axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

svi.minority.w <- svi.state_West %>%
  ggplot(aes(x = reorder(STATE, minority.rate), y = minority.rate, color = STATE, fill = STATE)) + 
  geom_bar(stat = "identity") + 
  coord_flip()+
  ylim(0, 0.8) +
  xlab("State") +
  ylab("Minority Rate") + 
  ggtitle("Minority Rate in West US") + 
  theme_minimal() +
  theme(legend.position = 'none', # remove legend,
        axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

Plotting disabled rate

svi.disabled.mw <- svi.state_Midwest %>%
  ggplot(aes(x = reorder(STATE, disabled.rate), y = disabled.rate, color = STATE, fill = STATE)) + 
  geom_bar(stat = "identity") + 
  coord_flip()+
  ylim(0, 0.2) +
  xlab("State") +
  ylab("Disabled Rate") + 
  ggtitle("Disabled Rate in Midwest US") + 
  theme_minimal() +
  theme(legend.position = 'none', # remove legend,
        axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

svi.disabled.ne <- svi.state_Northeast %>%
  ggplot(aes(x = reorder(STATE, disabled.rate), y = disabled.rate, color = STATE, fill = STATE)) + 
  geom_bar(stat = "identity") + 
  coord_flip()+
  ylim(0, 0.2) +
  xlab("State") +
  ylab("Disabled Rate") + 
  ggtitle("Disabled Rate in Northeast US") + 
  theme_minimal() +
  theme(legend.position = 'none', # remove legend,
        axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

svi.disabled.s <- svi.state_South %>%
  ggplot(aes(x = reorder(STATE, disabled.rate), y = disabled.rate, color = STATE, fill = STATE)) + 
  geom_bar(stat = "identity") + 
  coord_flip()+
  ylim(0, 0.2) +
  xlab("State") +
  ylab("Disabled Rate") + 
  ggtitle("Disabled Rate in South US") + 
  theme_minimal() +
  theme(legend.position = 'none', # remove legend,
        axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

svi.disabled.w <- svi.state_West %>%
  ggplot(aes(x = reorder(STATE, disabled.rate), y = disabled.rate, color = STATE, fill = STATE)) + 
  geom_bar(stat = "identity") + 
  coord_flip()+
  ylim(0, 0.2) +
  xlab("State") +
  ylab("Disabled Rate") + 
  ggtitle("Disabled Rate in West US") + 
  theme_minimal() +
  theme(legend.position = 'none', # remove legend,
        axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

Plotting elderly rate

svi.old.mw <- svi.state_Midwest %>%
  ggplot(aes(x = reorder(STATE, age65.rate), y = age65.rate, color = STATE, fill = STATE)) + 
  geom_bar(stat = "identity") + 
  coord_flip()+
  ylim(0, 0.2) +
  xlab("State") +
  ylab("Age over 65 Rate") + 
  ggtitle("Age over 65 Rate in Midwest US") + 
  theme_minimal() +
  theme(legend.position = 'none', # remove legend,
        axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

svi.old.ne <- svi.state_Northeast %>%
  ggplot(aes(x = reorder(STATE, age65.rate), y = age65.rate, color = STATE, fill = STATE)) + 
  geom_bar(stat = "identity") + 
  coord_flip()+
  ylim(0, 0.2) +
  xlab("State") +
  ylab("Age over 65 Rate") + 
  ggtitle("Age over 65 Rate in Northeast US") + 
  theme_minimal() +
  theme(legend.position = 'none', # remove legend,
        axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

svi.old.s <- svi.state_South %>%
  ggplot(aes(x = reorder(STATE, age65.rate), y = age65.rate, color = STATE, fill = STATE)) + 
  geom_bar(stat = "identity") + 
  coord_flip()+
  ylim(0, 0.2) +
  xlab("State") +
  ylab("Age over 65 Rate") + 
  ggtitle("Age over 65 Rate in South US") + 
  theme_minimal() +
  theme(legend.position = 'none', # remove legend,
        axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

svi.old.w <- svi.state_West %>%
  ggplot(aes(x = reorder(STATE, age65.rate), y = age65.rate, color = STATE, fill = STATE)) + 
  geom_bar(stat = "identity") + 
  coord_flip()+
  ylim(0, 0.2) +
  xlab("State") +
  ylab("Age over 65 Rate") + 
  ggtitle("Age over 65 Rate in West US") + 
  theme_minimal() +
  theme(legend.position = 'none', # remove legend,
        axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

Plotting per capita income

svi.pci.mw <- svi.state_Midwest %>%
  ggplot(aes(x = reorder(STATE, per.capita), y = per.capita, color = STATE, fill = STATE)) + 
  geom_bar(stat = "identity") + 
  coord_flip()+
  ylim(0, 55000) +
  xlab("State") +
  ylab("Per Capita Income") + 
  ggtitle("Per Capita Income in Midwest US") + 
  theme_minimal() +
  theme(legend.position = 'none', # remove legend,
        axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

svi.pci.ne <- svi.state_Northeast %>%
  ggplot(aes(x = reorder(STATE, per.capita), y = per.capita, color = STATE, fill = STATE)) + 
  geom_bar(stat = "identity") + 
  coord_flip()+
  ylim(0, 55000) +
  xlab("State") +
  ylab("Per Capita Income") + 
  ggtitle("Per Capita Income in Northeast US") + 
  theme_minimal() +
  theme(legend.position = 'none', # remove legend,
        axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

svi.pci.s <- svi.state_South %>%
  ggplot(aes(x = reorder(STATE, per.capita), y = per.capita, color = STATE, fill = STATE)) + 
  geom_bar(stat = "identity") + 
  coord_flip()+
  ylim(0, 55000) +
  xlab("State") +
  ylab("Per Capita Income") + 
  ggtitle("Per Capita Income in South US") + 
  theme_minimal() +
  theme(legend.position = 'none', # remove legend,
        axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style

svi.pci.w <- svi.state_West %>%
  ggplot(aes(x = reorder(STATE, per.capita), y = per.capita, color = STATE, fill = STATE)) + 
  geom_bar(stat = "identity") + 
  coord_flip()+
  geom_text(x=13, y =10800, label = "New Mexico data missing",
            size = 5, color = "black")+
  ylim(0, 55000) +
  xlab("State") +
  ylab("Per Capita Income") + 
  ggtitle("Per Capita Income in West US") + 
  theme_minimal() +
  theme(legend.position = 'none', # remove legend,
        axis.title = element_text(size=12), # change axis style
        legend.title = element_text(size=12), # change legend style
        plot.title = element_text(size=14)) # change title style
### New Mexico data missing

linear regression analysis

Create SVI & mortality rates data frame with mortality rates calculated by cumulative death over cumulative confirmed

svi.analysis <- svi.state %>%
  select(c(STATE, poverty.rate, per.capita, no.HS.rate, age65.rate, disabled.rate, minority.rate, uninsured.rate)) %>%
  mutate(state = STATE)

covid.mortality <- covid.time.series %>%
  filter(date == as.Date("2020-10-24")) %>%
  mutate(mortality.avg = (cumulative.death.cases / cumulative.confirmed.cases)*100) %>%
  select(state, mortality.avg)

svi.mortality <- svi.analysis %>%
  inner_join(covid.mortality, by = "state") %>%
  select(-state)

# log transformation of per capita income 
svi.mortality <- svi.mortality %>%
  mutate(pci.log = log10(per.capita))

Create mobility & active rates data frame with active rates calculated by active cases over total population with a lag of 5 days

library(lubridate)

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

# making changing rates of covid active data 
covid.active.diff <- covid.time.series %>% 
  select(c(state, date, active.rates)) %>%
  group_by(state) %>%
  mutate(active.changes = active.rates - lag(active.rates))

# potential complications of different baselines for each day in the week
# making changing rates of mobility data 
mobility.change <- StateMonthMobility %>%
  mutate(check.parks = parks_percent_change_from_baseline, 
         check.transit = transit_stations_percent_change_from_baseline) %>% 
#  fill(parks_percent_change_from_baseline) %>%
#  fill(transit_stations_percent_change_from_baseline) %>%
  mutate(retail.recration.diff = retail_and_recreation_percent_change_from_baseline - lag(retail_and_recreation_percent_change_from_baseline),
         grocery.pharmacy.diff = grocery_and_pharmacy_percent_change_from_baseline - lag(grocery_and_pharmacy_percent_change_from_baseline),
         parks.diff = parks_percent_change_from_baseline - lag(parks_percent_change_from_baseline),
         transit.diff = transit_stations_percent_change_from_baseline - lag(transit_stations_percent_change_from_baseline),
         workplaces.diff = workplaces_percent_change_from_baseline - lag(workplaces_percent_change_from_baseline),
         residential.diff = residential_percent_change_from_baseline - lag(residential_percent_change_from_baseline)) %>% 
  mutate(parks.diff = ifelse(is.na(check.parks), NA, parks.diff),
         transit.diff = ifelse(is.na(check.transit), NA, transit.diff))


sum(is.na(mobility.change$parks.diff))

## [1] 530

sum(is.na(mobility.change$parks_percent_change_from_baseline))

## [1] 440

# parks 440 NA
# transit 135 NA


mobility.active <- merge(mobility.change, covid.active.diff, by.x = c("sub_region_1", "date"), by.y = c("state", "date"))
# active rates df
mobility.active <- mobility.active %>%
  mutate(state = sub_region_1) %>%
  select(c(state, date, retail_and_recreation_percent_change_from_baseline, grocery_and_pharmacy_percent_change_from_baseline, parks_percent_change_from_baseline, transit_stations_percent_change_from_baseline, workplaces_percent_change_from_baseline, residential_percent_change_from_baseline, active.rates)) %>% 
  filter(date >= as.Date("2020-04-13")) %>% 
  mutate(month = months(date), 
         week = format(date, format = '%W'),
         # the following code converts to biweekly
         biweek = 14 * (as.numeric(date - min(date)) %/% 14) + min(date)) %>%
  group_by(state, biweek) %>%
  summarise(active.avg = mean(active.rates),
            retail.recreation.avg = mean(retail_and_recreation_percent_change_from_baseline),
            grocery.pharmacy.avg = mean(grocery_and_pharmacy_percent_change_from_baseline), 
            parks.avg = mean(parks_percent_change_from_baseline, na.rm=TRUE), 
            transit.avg = mean(transit_stations_percent_change_from_baseline, na.rm=TRUE),
            workplace.avg = mean(workplaces_percent_change_from_baseline),
            residential.avg = mean(residential_percent_change_from_baseline)) 

## `summarise()` regrouping output by 'state' (override with `.groups` argument)

# need to marking off public holidays 
# change to bi-weekly later
# use SVI to pick out top n and bottom n states 
# within these states, do state by state mobility ~ covid lm 

Linear regression

svi.crit = 10

# poverty rates
svi.pov.top <- svi.mortality %>%
  slice_max(poverty.rate, n = svi.crit)
svi.pov.bottom <- svi.mortality %>%
  slice_min(poverty.rate, n = svi.crit)

svi.pov.top.states <- svi.pov.top$STATE
svi.pov.bottom.states <- svi.pov.bottom$STATE 

pov.mortality.lm <- lm(mortality.avg ~ poverty.rate, svi.mortality)
pov.mortality.sum <- summary(pov.mortality.lm)
pov.mortality.top.lm <- lm(mortality.avg ~ poverty.rate, svi.pov.top)
pov.mortality.top.sum <- summary(pov.mortality.top.lm)
pov.mortality.bottom.lm <- lm(mortality.avg ~ poverty.rate, svi.pov.bottom)
pov.mortality.bottom.sum <- summary(pov.mortality.bottom.lm)


# per capita 
svi.pci.top <- svi.mortality %>%
  slice_max(per.capita, n = svi.crit)
svi.pci.bottom <- svi.mortality %>%
  slice_min(per.capita, n = svi.crit)

svi.pci.top.states <- svi.pci.top$STATE
svi.pci.bottom.states <- svi.pci.bottom$STATE 

pci.log.mortality.lm <- lm(mortality.avg ~ pci.log, svi.mortality)
pci.log.mortality.sum <- summary(pci.log.mortality.lm)
pci.log.mortality.top.lm <- lm(mortality.avg ~ pci.log, svi.pci.top)
pci.log.mortality.top.sum <- summary(pci.log.mortality.top.lm)
pci.log.mortality.bottom.lm <- lm(mortality.avg ~ pci.log, svi.pci.bottom)
pci.log.mortality.bottom.sum <- summary(pci.log.mortality.bottom.lm)
# pci is significant for all states but not top and bottom states 


# no HS rates 
svi.noHS.top <- svi.mortality %>%
  slice_max(no.HS.rate, n = svi.crit)
svi.noHS.bottom <- svi.mortality %>%
  slice_min(no.HS.rate, n = svi.crit)

svi.noHS.top.states <- svi.noHS.top$STATE
svi.noHS.bottom.states <- svi.noHS.bottom$STATE 



noHS.mortality.top.lm <- lm(mortality.avg ~ no.HS.rate, svi.noHS.top)
noHS.mortality.top.sum <- summary(noHS.mortality.top.lm)
noHS.mortality.bottom.lm <- lm(mortality.avg ~ no.HS.rate, svi.noHS.bottom)
noHS.mortality.bottom.sum <- summary(noHS.mortality.bottom.lm)



# age 65/elderly rates 
svi.age.top <- svi.mortality %>%
  slice_max(age65.rate, n = svi.crit)
svi.age.bottom <- svi.mortality %>%
  slice_min(age65.rate, n = svi.crit)

svi.age.top.states <- svi.age.top$STATE
svi.age.bottom.states <- svi.age.bottom$STATE

age.mortality.top.lm <- lm(mortality.avg ~ age65.rate, svi.age.top)
age.mortality.top.sum <- summary(age.mortality.top.lm)
age.mortality.bottom.lm <- lm(mortality.avg ~ age65.rate, svi.age.bottom)
age.mortality.bottom.sum <- summary(age.mortality.bottom.lm)
# age bottom


# disabled rate
svi.disabled.top <- svi.mortality %>%
  slice_max(disabled.rate, n = svi.crit)
svi.disabled.bottom <- svi.mortality %>%
  slice_min(disabled.rate, n = svi.crit)

svi.disabled.top.states <- svi.disabled.top$STATE
svi.disabled.bottom.states <- svi.disabled.bottom$STATE 

disabled.mortality.top.lm <- lm(mortality.avg ~ disabled.rate, svi.disabled.top)
disabled.mortality.top.sum <- summary(disabled.mortality.top.lm)
disabled.mortality.bottom.lm <- lm(mortality.avg ~ disabled.rate, svi.disabled.bottom)
disabled.mortality.bottom.sum <- summary(disabled.mortality.bottom.lm)


# minority 
svi.minority.top <- svi.mortality %>%
        slice_max(minority.rate, n = svi.crit)
svi.minority.bottom <- svi.mortality %>%
        slice_min(minority.rate, n = svi.crit)

svi.minority.top.states <- svi.minority.top$STATE
svi.minority.bottom.states <- svi.minority.bottom$STATE

minority.mortality.top.lm <- lm(mortality.avg ~ minority.rate, svi.minority.top)
minority.mortality.top.sum <- summary(minority.mortality.top.lm)
minority.mortality.bottom.lm <- lm(mortality.avg ~ minority.rate, svi.minority.bottom)
minority.mortality.bottom.sum <- summary(minority.mortality.bottom.lm)
# minority bottom


# uninsured rates 
svi.uninsured.top <- svi.mortality %>%
        slice_max(uninsured.rate, n = svi.crit)
svi.uninsured.bottom <- svi.mortality %>%
        slice_min(uninsured.rate, n = svi.crit)

svi.uninsured.top.states <- svi.uninsured.top$STATE
svi.uninsured.bottom.states <- svi.uninsured.bottom$STATE

uninsured.mortality.lm <- lm(mortality.avg ~ uninsured.rate, svi.mortality)
uninsured.mortality.sum <- summary(uninsured.mortality.lm)
uninsured.mortality.top.lm <- lm(mortality.avg ~ uninsured.rate, svi.uninsured.top)
uninsured.mortality.top.sum <- summary(uninsured.mortality.top.lm)
uninsured.mortality.bottom.lm <- lm(mortality.avg ~ uninsured.rate, svi.uninsured.bottom)
uninsured.mortality.bottom.sum <- summary(uninsured.mortality.bottom.lm)
# uninsured also significant 

# pci & uninsured
pci.log.uninsured.mortality.lm <- lm(mortality.avg ~ pci.log + uninsured.rate, svi.mortality)
pci.log.uninsured.mortality.sum <- summary(pci.log.uninsured.mortality.lm)
pci.uninsured.interaction.mortality.lm <- lm(mortality.avg ~ pci.log * uninsured.rate, svi.mortality)
pci.uninsured.interaction.mortality.sum <- summary(pci.uninsured.interaction.mortality.lm)
pci.uninsured <- svi.mortality %>%
  filter(!is.na(pci.log))
cor(pci.uninsured$pci.log, pci.uninsured$uninsured.rate)

## [1] -0.4745319

# cor < 0.7 

svi.top.states <- table(c(svi.pov.top.states, svi.pci.bottom.states, svi.noHS.top.states, svi.age.top.states, svi.disabled.top.states, svi.minority.top.states, svi.uninsured.top.states)) 

svi.top.overlap.states <- sort(svi.top.states, decreasing=T)[1:5]

svi.bottom.states <- table(c(svi.pov.bottom.states, svi.pci.top.states, svi.noHS.bottom.states, svi.age.bottom.states, svi.disabled.bottom.states, svi.minority.bottom.states, svi.uninsured.bottom.states)) 
svi.bottom.overlap.states <- sort(svi.bottom.states, decreasing=T)[1:6]

# top and bottom states for pci and uninsured 
svi.sigl.crit = 5
svi.pci.sig.top <- svi.mortality %>%
  slice_max(per.capita, n = svi.sigl.crit)
svi.pci.sig.top.states = svi.pci.sig.top$STATE
#  *"Connecticut"*   "Massachusetts" "New Jersey"    *"Maryland"*      "New Hampshire"
svi.pci.sig.bottom <- svi.mortality %>%
  slice_min(per.capita, n = svi.sigl.crit)
svi.pci.sig.bottom.states = svi.pci.sig.bottom$STATE
# *"Mississippi"*   *"West Virginia"* "Arkansas"      "Idaho"         *"Alabama"* 

svi.sig.uninsured.top <- svi.mortality %>%
        slice_max(uninsured.rate, n = svi.sigl.crit)
svi.sig.uninsured.top.states = svi.sig.uninsured.top$STATE
# "Texas"    "Alaska"   "Oklahoma" "Georgia"  "Florida" 
svi.sig.uninsured.bottom <- svi.mortality %>%
        slice_min(uninsured.rate, n = svi.sigl.crit)
svi.sig.uninsured.bottom.states = svi.sig.uninsured.bottom$STATE
# "Massachusetts" "Hawaii"        "Vermont"       *"Minnesota"*     "Iowa"   

svi states to covid ~ active rates

Top Five States

Mississippi

#mobility.active.svi.top <- mobility.active %>%
#  filter(state %in% c("Mississippi", "West Virginia", "Alabama", "Kentucky", "Louisiana"))

mobility.active.Mississippi <- mobility.active %>%
  filter(state == "Mississippi")

retail.active.Mississippi.lm <- lm(active.avg ~ retail.recreation.avg, mobility.active.Mississippi)
retail.active.Mississippi.sum <- summary(retail.active.Mississippi.lm)

grocery.active.Mississippi.lm <- lm(active.avg ~ grocery.pharmacy.avg, mobility.active.Mississippi)
grocery.active.Mississippi.sum <- summary(grocery.active.Mississippi.lm)

parks.active.Mississippi.lm <- lm(active.avg ~ parks.avg, mobility.active.Mississippi)
parks.active.Mississippi.sum <- summary(parks.active.Mississippi.lm)

transit.active.Mississippi.lm <- lm(active.avg ~ transit.avg, mobility.active.Mississippi)
transit.active.Mississippi.sum <- summary(transit.active.Mississippi.lm)

workplace.active.Mississippi.lm <- lm(active.avg ~ workplace.avg, mobility.active.Mississippi)
workplace.active.Mississippi.sum <- summary(workplace.active.Mississippi.lm)

residential.active.Mississippi.lm <- lm(active.avg ~ residential.avg, mobility.active.Mississippi)
residential.active.Mississippi.sum <- summary(residential.active.Mississippi.lm)

# significant results: none :(

West Virginia

mobility.active.West.Virginia <- mobility.active %>%
  filter(state == "West Virginia")

retail.active.West.Virginia.lm <- lm(active.avg ~ retail.recreation.avg, mobility.active.West.Virginia)
retail.active.West.Virginia.sum <- summary(retail.active.West.Virginia.lm)

grocery.active.West.Virginia.lm <- lm(active.avg ~ grocery.pharmacy.avg, mobility.active.West.Virginia)
grocery.active.West.Virginia.sum <- summary(grocery.active.West.Virginia.lm)

parks.active.West.Virginia.lm <- lm(active.avg ~ parks.avg, mobility.active.West.Virginia)
parks.active.West.Virginia.sum <- summary(parks.active.West.Virginia.lm)

transit.active.West.Virginia.lm <- lm(active.avg ~ transit.avg, mobility.active.West.Virginia)
transit.active.West.Virginia.sum <- summary(transit.active.West.Virginia.lm)

workplace.active.West.Virginia.lm <- lm(active.avg ~ workplace.avg, mobility.active.West.Virginia)
workplace.active.West.Virginia.sum <- summary(workplace.active.West.Virginia.lm)

residential.active.West.Virginia.lm <- lm(active.avg ~ residential.avg, mobility.active.West.Virginia)
residential.active.West.Virginia.sum <- summary(residential.active.West.Virginia.lm)

# significant results: workplace.active.West.Virginia.sum(**), residential.active.West.Virginia.sum(*)

Alabama

mobility.active.Alabama <- mobility.active %>%
  filter(state == "Alabama")

retail.active.Alabama.lm <- lm(active.avg ~ retail.recreation.avg, mobility.active.Alabama)
retail.active.Alabama.sum <- summary(retail.active.Alabama.lm)

grocery.active.Alabama.lm <- lm(active.avg ~ grocery.pharmacy.avg, mobility.active.Alabama)
grocery.active.Alabama.sum <- summary(grocery.active.Alabama.lm)

parks.active.Alabama.lm <- lm(active.avg ~ parks.avg, mobility.active.Alabama)
parks.active.Alabama.sum <- summary(parks.active.Alabama.lm)

transit.active.Alabama.lm <- lm(active.avg ~ transit.avg, mobility.active.Alabama)
transit.active.Alabama.sum <- summary(transit.active.Alabama.lm)

workplace.active.Alabama.lm <- lm(active.avg ~ workplace.avg, mobility.active.Alabama)
workplace.active.Alabama.sum <- summary(workplace.active.Alabama.lm)

residential.active.Alabama.lm <- lm(active.avg ~ residential.avg, mobility.active.Alabama)
residential.active.Alabama.sum <- summary(residential.active.Alabama.lm)

# significant results: grocery.active.Alabama.sum(*), workplace.active.Alabama.sum(*), residential.active.Alabama.sum (***)

Kentucky

mobility.active.Kentucky  <- mobility.active %>%
  filter(state == "Kentucky")

retail.active.Kentucky.lm <- lm(active.avg ~ retail.recreation.avg, mobility.active.Kentucky)
retail.active.Kentucky.sum <- summary(retail.active.Kentucky.lm)

grocery.active.Kentucky.lm <- lm(active.avg ~ grocery.pharmacy.avg, mobility.active.Kentucky)
grocery.active.Kentucky.sum <- summary(grocery.active.Kentucky.lm)

parks.active.Kentucky.lm <- lm(active.avg ~ parks.avg, mobility.active.Kentucky)
parks.active.Kentucky.sum <- summary(parks.active.Kentucky.lm)

transit.active.Kentucky.lm <- lm(active.avg ~ transit.avg, mobility.active.Kentucky)
transit.active.Kentucky.sum <- summary(transit.active.Kentucky.lm)

workplace.active.Kentucky.lm <- lm(active.avg ~ workplace.avg, mobility.active.Kentucky)
workplace.active.Kentucky.sum <- summary(workplace.active.Kentucky.lm)

residential.active.Kentucky.lm <- lm(active.avg ~ residential.avg, mobility.active.Kentucky)
residential.active.Kentucky.sum <- summary(residential.active.Kentucky.lm)

# significant results: transit.active.Kentucky.sum(*), workplace.active.Kentucky.sum(**), residential.active.Kentucky.sum(*)

Louisiana

mobility.active.Louisiana  <- mobility.active %>%
  filter(state == "Louisiana")

retail.active.Louisiana.lm <- lm(active.avg ~ retail.recreation.avg, mobility.active.Louisiana)
retail.active.Louisiana.sum <- summary(retail.active.Louisiana.lm)

grocery.active.Louisiana.lm <- lm(active.avg ~ grocery.pharmacy.avg, mobility.active.Louisiana)
grocery.active.Louisiana.sum <- summary(grocery.active.Louisiana.lm)

parks.active.Louisiana.lm <- lm(active.avg ~ parks.avg, mobility.active.Louisiana)
parks.active.Louisiana.sum <- summary(parks.active.Louisiana.lm)

transit.active.Louisiana.lm <- lm(active.avg ~ transit.avg, mobility.active.Louisiana)
transit.active.Louisiana.sum <- summary(transit.active.Louisiana.lm)

workplace.active.Louisiana.lm <- lm(active.avg ~ workplace.avg, mobility.active.Louisiana)
workplace.active.Louisiana.sum <- summary(workplace.active.Louisiana.lm)

residential.active.Louisiana.lm <- lm(active.avg ~ residential.avg, mobility.active.Louisiana)
residential.active.Louisiana.sum <- summary(residential.active.Louisiana.lm)

# significant results: none :(

Bottom Six States

Minnesota

mobility.active.Minnesota <- mobility.active %>%
  filter(state == "Minnesota")

#retail and recreation 
retail.active.Minnesota.lm <- lm(active.avg ~ retail.recreation.avg, mobility.active.Minnesota)
retail.active.Minnesota.sum <- summary(retail.active.Minnesota.lm)


#grocery 
grocery.active.Minnesota.lm <- lm(active.avg ~ grocery.pharmacy.avg, mobility.active.Minnesota)
grocery.active.Minnesota.sum <- summary(grocery.active.Minnesota.lm)

#parks
parks.active.Minnesota.lm <- lm(active.avg ~ parks.avg, mobility.active.Minnesota)
parks.active.Minnesota.sum <- summary(parks.active.Minnesota.lm)


#transit
transit.active.Minnesota.lm <- lm(active.avg ~ transit.avg, mobility.active.Minnesota)
transit.active.Minnesota.sum <- summary(transit.active.Minnesota.lm)

#workplace
workplace.active.Minnesota.lm <- lm(active.avg ~ workplace.avg, mobility.active.Minnesota)
workplace.active.Minnesota.sum <- summary(workplace.active.Minnesota.lm)

#residential 
residential.active.Minnesota.lm <- lm(active.avg ~ residential.avg, mobility.active.Minnesota)
residential.active.Minnesota.sum <- summary(residential.active.Minnesota.lm)


#none significant for minnesota 

North Dakota

mobility.active.ND <- mobility.active %>%
  filter(state == "North Dakota")

#retail and recreation 
retail.active.ND.lm <- lm(active.avg ~ retail.recreation.avg, mobility.active.ND)
retail.active.ND.sum <- summary(retail.active.ND.lm)


#grocery 
grocery.active.ND.lm <- lm(active.avg ~ grocery.pharmacy.avg, mobility.active.ND)
grocery.active.ND.sum <- summary(grocery.active.ND.lm)

#parks
parks.active.ND.lm <- lm(active.avg ~ parks.avg, mobility.active.ND)
parks.active.ND.sum <- summary(parks.active.ND.lm)


#transit
transit.active.ND.lm <- lm(active.avg ~ transit.avg, mobility.active.ND)
transit.active.ND.sum <- summary(transit.active.ND.lm)

#workplace
workplace.active.ND.lm <- lm(active.avg ~ workplace.avg, mobility.active.ND)
workplace.active.ND.sum <- summary(workplace.active.ND.lm)

#residential 
residential.active.ND.lm <- lm(active.avg ~ residential.avg, mobility.active.ND)
residential.active.ND.sum <- summary(residential.active.ND.lm)

#grocery.pharmacy.avg significant (at a = 0.05) with p = 0.013257

Connecticut

mobility.active.Connecticut <- mobility.active %>%
  filter(state == "Connecticut")

retail.active.Connecticut.lm <- lm(active.avg ~ retail.recreation.avg, mobility.active.Connecticut)
retail.active.Connecticut.sum <- summary(retail.active.Connecticut.lm)

grocery.active.Connecticut.lm <- lm(active.avg ~ grocery.pharmacy.avg, mobility.active.Connecticut)
grocery.active.Connecticut.sum <- summary(grocery.active.Connecticut.lm)

parks.active.Connecticut.lm <- lm(active.avg ~ parks.avg, mobility.active.Connecticut)
parks.active.Connecticut.sum <- summary(parks.active.Connecticut.lm)

transit.active.Connecticut.lm <- lm(active.avg ~ transit.avg, mobility.active.Connecticut)
transit.active.Connecticut.sum <- summary(transit.active.Connecticut.lm)

workplace.active.Connecticut.lm <- lm(active.avg ~ workplace.avg, mobility.active.Connecticut)
workplace.active.Connecticut.sum <- summary(workplace.active.Connecticut.lm)

residential.active.Connecticut.lm <- lm(active.avg ~ residential.avg, mobility.active.Connecticut)
residential.active.Connecticut.sum <- summary(residential.active.Connecticut.lm)

#significant results: retail.active.Connecticut.sum (***), grocery.active.Connecticut.sum (*),transit.active.Connecticut.sum (***), workplace.active.Connecticut.sum (***), residential.active.Connecticut.sum (***)

Maryland

mobility.active.Maryland <- mobility.active %>%
  filter(state == "Maryland")

retail.active.Maryland.lm <- lm(active.avg ~ retail.recreation.avg, mobility.active.Maryland)
retail.active.Maryland.sum <- summary(retail.active.Maryland.lm)

grocery.active.Maryland.lm <- lm(active.avg ~ grocery.pharmacy.avg, mobility.active.Maryland)
grocery.active.Maryland.sum <- summary(grocery.active.Maryland.lm)

parks.active.Maryland.lm <- lm(active.avg ~ parks.avg, mobility.active.Maryland)
parks.active.Maryland.sum <- summary(parks.active.Maryland.lm)

transit.active.Maryland.lm <- lm(active.avg ~ transit.avg, mobility.active.Maryland)
transit.active.Maryland.sum <- summary(transit.active.Maryland.lm)

workplace.active.Maryland.lm <- lm(active.avg ~ workplace.avg, mobility.active.Maryland)
workplace.active.Maryland.sum <- summary(workplace.active.Maryland.lm)

residential.active.Maryland.lm <- lm(active.avg ~ residential.avg, mobility.active.Maryland)
residential.active.Maryland.sum <- summary(residential.active.Maryland.lm)

#significant results: retail.active.Maryland.sum (***), transit.active.Maryland.sum (*), workplace.active.Maryland.sum (***), workplace.active.Maryland.sum (***). 

New Hampshire

mobility.active.New.Hampshire <- mobility.active %>%
  filter(state == "New Hampshire")

retail.active.New.Hampshire.lm <- lm(active.avg ~ retail.recreation.avg, mobility.active.New.Hampshire)
retail.active.New.Hampshire.sum <- summary(retail.active.New.Hampshire.lm)


grocery.active.New.Hampshire.lm <- lm(active.avg ~ grocery.pharmacy.avg, mobility.active.New.Hampshire)
grocery.active.New.Hampshire.sum <- summary(grocery.active.New.Hampshire.lm)

parks.active.New.Hampshire.lm <- lm(active.avg ~ parks.avg, mobility.active.New.Hampshire)
parks.active.New.Hampshire.sum <- summary(parks.active.New.Hampshire.lm)

transit.active.New.Hampshire.lm <- lm(active.avg ~ transit.avg, mobility.active.New.Hampshire)
transit.active.New.Hampshire.sum <- summary(transit.active.New.Hampshire.lm)



workplace.active.New.Hampshire.lm <- lm(active.avg ~ workplace.avg, mobility.active.New.Hampshire)
workplace.active.New.Hampshire.sum <- summary(workplace.active.New.Hampshire.lm)



residential.active.New.Hampshire.lm <- lm(active.avg ~ residential.avg, mobility.active.New.Hampshire)
residential.active.New.Hampshire.sum <- summary(residential.active.New.Hampshire.lm)



# significant results:retail.active.New.Hampshire.sum(**),grocery.active.New.Hampshire.sum(**),parks.active.New.Hampshire.sum(*),
#transit.active.New.Hampshire.sum(**),workplace.active.New.Hampshire.sum(**),residential.active.New.Hampshire.sum(***)

Utah

mobility.active.Utah <- mobility.active %>%
  filter(state == "Utah")

retail.active.Utah.lm <- lm(active.avg ~ retail.recreation.avg, mobility.active.Utah)
retail.active.Utah.sum <- summary(retail.active.Utah.lm)



grocery.active.Utah.lm <- lm(active.avg ~ grocery.pharmacy.avg, mobility.active.Utah)
grocery.active.Utah.sum <- summary(grocery.active.Utah.lm)


parks.active.Utah.lm <- lm(active.avg ~ parks.avg, mobility.active.Utah)
parks.active.Utah.sum <- summary(parks.active.Utah.lm)



transit.active.Utah.lm <- lm(active.avg ~ transit.avg, mobility.active.Utah)
transit.active.Utah.sum <- summary(transit.active.Utah.lm)


workplace.active.Utah.lm <- lm(active.avg ~ workplace.avg, mobility.active.Utah)
workplace.active.Utah.sum <- summary(workplace.active.Utah.lm)


residential.active.Utah.lm <- lm(active.avg ~ residential.avg, mobility.active.Utah)
residential.active.Utah.sum <- summary(residential.active.Utah.lm)


# significant results:  transit.active.Utah.sum (*),workplace.active.Utah.sum(*),residental.active.Utah.sum(*)

Massachusetts

mobility.active.Massachusetts <- mobility.active %>%
  filter(state == "Massachusetts")

retail.active.Massachusetts.lm <- lm(active.avg ~ retail.recreation.avg, mobility.active.Massachusetts)
retail.active.Massachusetts.sum <- summary(retail.active.Massachusetts.lm)



grocery.active.Massachusetts.lm <- lm(active.avg ~ grocery.pharmacy.avg, mobility.active.Massachusetts)
grocery.active.Massachusetts.sum <- summary(grocery.active.Massachusetts.lm)


parks.active.Massachusetts.lm <- lm(active.avg ~ parks.avg, mobility.active.Massachusetts)
parks.active.Massachusetts.sum <- summary(parks.active.Massachusetts.lm)



transit.active.Massachusetts.lm <- lm(active.avg ~ transit.avg, mobility.active.Massachusetts)
transit.active.Massachusetts.sum <- summary(transit.active.Massachusetts.lm)


workplace.active.Massachusetts.lm <- lm(active.avg ~ workplace.avg, mobility.active.Massachusetts)
workplace.active.Massachusetts.sum <- summary(workplace.active.Massachusetts.lm)


residential.active.Massachusetts.lm <- lm(active.avg ~ residential.avg, mobility.active.Massachusetts)
residential.active.Massachusetts.sum <- summary(residential.active.Massachusetts.lm)


# significant results:  none

New Jersey

mobility.active.New.Jersey <- mobility.active %>%
  filter(state == "New Jersey")

retail.active.New.Jersey.lm <- lm(active.avg ~ retail.recreation.avg, mobility.active.New.Jersey)
retail.active.New.Jersey.sum <- summary(retail.active.New.Jersey.lm)



grocery.active.New.Jersey.lm <- lm(active.avg ~ grocery.pharmacy.avg, mobility.active.New.Jersey)
grocery.active.New.Jersey.sum <- summary(grocery.active.New.Jersey.lm)


parks.active.New.Jersey.lm <- lm(active.avg ~ parks.avg, mobility.active.New.Jersey)
parks.active.New.Jersey.sum <- summary(parks.active.New.Jersey.lm)



transit.active.New.Jersey.lm <- lm(active.avg ~ transit.avg, mobility.active.New.Jersey)
transit.active.New.Jersey.sum <- summary(transit.active.New.Jersey.lm)


workplace.active.New.Jersey.lm <- lm(active.avg ~ workplace.avg, mobility.active.New.Jersey)
workplace.active.New.Jersey.sum <- summary(workplace.active.New.Jersey.lm)


residential.active.New.Jersey.lm <- lm(active.avg ~ residential.avg, mobility.active.New.Jersey)
residential.active.New.Jersey.sum <- summary(residential.active.New.Jersey.lm)


# significant results:  retail.active.New.Jersey.sum(***), transit.active.New.Jersey.sum(***), workplace.active.New.Jersey.sum (***), residential.active.New.Jersey.sum(***)

Texas

mobility.active.Texas <- mobility.active %>%
        filter(state == "Texas")
#retail and recreation 
retail.active.Texas.lm <- lm(active.avg ~ retail.recreation.avg, mobility.active.Texas)
retail.active.Texas.sum <- summary(retail.active.Texas.lm)
#grocery 
grocery.active.Texas.lm <- lm(active.avg ~ grocery.pharmacy.avg, mobility.active.Texas)
grocery.active.Texas.sum <- summary(grocery.active.Texas.lm)
#parks
parks.active.Texas.lm <- lm(active.avg ~ parks.avg, mobility.active.Texas)
parks.active.Texas.sum <- summary(parks.active.Texas.lm)
#transit
transit.active.Texas.lm <- lm(active.avg ~ transit.avg, mobility.active.Texas)
transit.active.Texas.sum <- summary(transit.active.Texas.lm)
#workplace
workplace.active.Texas.lm <- lm(active.avg ~ workplace.avg, mobility.active.Texas)
workplace.active.Texasa.sum <- summary(workplace.active.Texas.lm)
#residential 
residential.active.Texas.lm <- lm(active.avg ~ residential.avg, mobility.active.Texas)
residential.active.Texas.sum <- summary(residential.active.Texas.lm)

# significant results:  transit.active.Texas.sum(*)

Arkansas

mobility.active.Arkansas <- mobility.active %>%        
  filter(state == "Arkansas")
#retail and recreation 
retail.active.Arkansas.lm <- lm(active.avg ~ retail.recreation.avg, mobility.active.Arkansas)
retail.active.Arkansas.sum <- summary(retail.active.Arkansas.lm)
#grocery 
grocery.active.Arkansas.lm <- lm(active.avg ~ grocery.pharmacy.avg, mobility.active.Arkansas)
grocery.active.Arkansas.sum <- summary(grocery.active.Arkansas.lm)
#parks 
parks.active.Arkansas.lm <- lm(active.avg ~ parks.avg, mobility.active.Arkansas) 
parks.active.Arkansas.sum <- summary(parks.active.Arkansas.lm)

#transit
transit.active.Arkansas.lm <- lm(active.avg ~ transit.avg, mobility.active.Arkansas)
transit.active.Arkansas.sum <- summary(transit.active.Arkansas.lm)
#workplace
workplace.active.Arkansas.lm <- lm(active.avg ~ workplace.avg, mobility.active.Arkansas)
workplace.active.Arkansas.sum <- summary(workplace.active.Arkansas.lm)

#residential 
residential.active.Arkansas.lm <- lm(active.avg ~ residential.avg, mobility.active.Arkansas)
residential.active.Arkansas.sum <- summary(residential.active.Arkansas.lm)

# significant results: transit.active.Arkansas.sum(**), workplace.active.Arkansas.sum(*), residential.active.Arkansas.sum(**)

Georgia

mobility.active.Georgia <- mobility.active %>%
  filter(state == "Georgia")

retail.active.Georgia.lm <- lm(active.avg ~ retail.recreation.avg, mobility.active.Georgia)
retail.active.Georgia.sum <- summary(retail.active.Georgia.lm)



grocery.active.Goergia.lm <- lm(active.avg ~ grocery.pharmacy.avg, mobility.active.Georgia)
grocery.active.Goergia.sum <- summary(grocery.active.Goergia.lm)


parks.active.Goergia.lm <- lm(active.avg ~ parks.avg, mobility.active.Georgia)
parks.active.Goergia.sum <- summary(parks.active.Goergia.lm)



transit.active.Goergia.lm <- lm(active.avg ~ transit.avg, mobility.active.Georgia)
transit.active.Goergia.sum <- summary(transit.active.Goergia.lm)


workplace.active.Goergia.lm <- lm(active.avg ~ workplace.avg, mobility.active.Georgia)
workplace.active.Goergia.sum <- summary(workplace.active.Goergia.lm)


residential.active.Goergia.lm <- lm(active.avg ~ residential.avg, mobility.active.Georgia)
residential.active.Goergia.sum <- summary(residential.active.Goergia.lm)

# significant results:  transit.active.Goergia.sum(**), workplace.active.Goergia.sum(**), residential.active.Goergia.sum(***)

Idaho

mobility.active.Idaho <- mobility.active %>%        
  filter(state == "Idaho")
#retail and recreation 
retail.active.Idaho.lm <- lm(active.avg ~ retail.recreation.avg, mobility.active.Idaho)
retail.active.Idaho.sum <- summary(retail.active.Idaho.lm)
#grocery 
grocery.active.Idaho.lm <- lm(active.avg ~ grocery.pharmacy.avg, mobility.active.Idaho)
grocery.active.Idaho.sum <- summary(grocery.active.Idaho.lm)
#parks
parks.active.Idaho.lm <- lm(active.avg ~ parks.avg, mobility.active.Idaho)
parks.active.Idaho.sum <- summary(parks.active.Idaho.lm)

#transit
transit.active.Idaho.lm <- lm(active.avg ~ transit.avg, mobility.active.Idaho)
transit.active.Idaho.sum <- summary(transit.active.Idaho.lm)
#workplace
workplace.active.Idaho.lm <- lm(active.avg ~ workplace.avg, mobility.active.Idaho)
workplace.active.Idaho.sum <- summary(workplace.active.Idaho.lm)

#residential 
residential.active.Idaho.lm <- lm(active.avg ~ residential.avg, mobility.active.Idaho)
residential.active.Idaho.sum <- summary(residential.active.Idaho.lm)


# significant results:  workplace.active.Idaho.sum(*), residential.active.Idaho.sum(*)

Florida

mobility.active.Florida <- mobility.active %>%
  filter(state == "Florida")

retail.active.Florida.lm <- lm(active.avg ~ retail.recreation.avg, mobility.active.Florida)
retail.active.Florida.sum <- summary(retail.active.Florida.lm)



grocery.active.Florida.lm <- lm(active.avg ~ grocery.pharmacy.avg, mobility.active.Florida)
grocery.active.Florida.sum <- summary(grocery.active.Florida.lm)


parks.active.Florida.lm <- lm(active.avg ~ parks.avg, mobility.active.Florida)
parks.active.Florida.sum <- summary(parks.active.Florida.lm)



transit.active.Florida.lm <- lm(active.avg ~ transit.avg, mobility.active.Florida)
transit.active.Florida.sum <- summary(transit.active.Florida.lm)


workplace.active.Florida.lm <- lm(active.avg ~ workplace.avg, mobility.active.Florida)
workplace.active.Florida.sum <- summary(workplace.active.Florida.lm)


residential.active.Florida.lm <- lm(active.avg ~ residential.avg, mobility.active.Florida)
residential.active.Florida.sum <- summary(residential.active.Florida.lm)

# significant results:  transit.active.Florida.sum(*), workplace.active.Florida.sum(*), residential.active.Florida.sum(***)

Oklahoma

mobility.active.Oklahoma <- mobility.active %>%
        filter(state == "Oklahoma")
#retail and recreation 
retail.active.Oklahoma.lm <- lm(active.avg ~ retail.recreation.avg, mobility.active.Oklahoma)
retail.active.Oklahoma.sum <- summary(retail.active.Oklahoma.lm)


#grocery 
grocery.active.Oklahoma.lm <- lm(active.avg ~ grocery.pharmacy.avg, mobility.active.Oklahoma)
grocery.active.Oklahoma.sum <- summary(grocery.active.Oklahoma.lm)
#parks
parks.active.Oklahoma.lm <- lm(active.avg ~ parks.avg, mobility.active.Oklahoma)
parks.active.Oklahoma.sum <- summary(parks.active.Oklahoma.lm)


#transit
transit.active.Oklahoma.lm <- lm(active.avg ~ transit.avg, mobility.active.Oklahoma)
transit.active.Oklahoma.sum <- summary(transit.active.Oklahoma.lm)
#workplace
workplace.active.Oklahoma.lm <- lm(active.avg ~ workplace.avg, mobility.active.Oklahoma)
workplace.active.Oklahoma.sum <- summary(workplace.active.Oklahoma.lm)

#residential 
residential.active.Oklahoma.lm <- lm(active.avg ~ residential.avg, mobility.active.Oklahoma)
residential.active.Oklahoma.sum <- summary(residential.active.Oklahoma.lm)

# significant results:  grocery.active.Oklahoma.sum(*), workplace.active.Oklahoma.sum(*), residential.active.Oklahoma.sum(*)

Hawaii

mobility.active.Hawaii <- mobility.active %>%
  filter(state == "Hawaii")

retail.active.Hawaii.lm <- lm(active.avg ~ retail.recreation.avg, mobility.active.Hawaii)
retail.active.Hawaii.sum <- summary(retail.active.Hawaii.lm)

grocery.active.Hawaii.lm <- lm(active.avg ~ grocery.pharmacy.avg, mobility.active.Hawaii)
grocery.active.Hawaii.sum <- summary(grocery.active.Hawaii.lm)

parks.active.Hawaii.lm <- lm(active.avg ~ parks.avg, mobility.active.Hawaii)
parks.active.Hawaii.sum <- summary(parks.active.Hawaii.lm)

transit.active.Hawaii.lm <- lm(active.avg ~ transit.avg, mobility.active.Hawaii)
transit.active.Hawaii.sum <- summary(transit.active.Hawaii.lm)

workplace.active.Hawaii.lm <- lm(active.avg ~ workplace.avg, mobility.active.Hawaii)
workplace.active.Hawaii.sum <- summary(workplace.active.Hawaii.lm)

residential.active.Hawaii.lm <- lm(active.avg ~ residential.avg, mobility.active.Hawaii)
residential.active.Hawaii.sum <- summary(residential.active.Hawaii.lm)
#significant results: grocery.active.Hawaii.sum (*),parks.active.Hawaii.sum (*)

Vermont

mobility.active.Vermont <- mobility.active %>%
  filter(state == "Vermont")

retail.active.Vermont.lm <- lm(active.avg ~ retail.recreation.avg, mobility.active.Vermont)
retail.active.Vermont.sum <- summary(retail.active.Vermont.lm)

grocery.active.Vermont.lm <- lm(active.avg ~ grocery.pharmacy.avg, mobility.active.Vermont)
grocery.active.Vermont.sum <- summary(grocery.active.Vermont.lm)

parks.active.Vermont.lm <- lm(active.avg ~ parks.avg, mobility.active.Vermont)
parks.active.Vermont.sum <- summary(parks.active.Vermont.lm)

transit.active.Vermont.lm <- lm(active.avg ~ transit.avg, mobility.active.Vermont)
transit.active.Vermont.sum <- summary(transit.active.Vermont.lm)

workplace.active.Vermont.lm <- lm(active.avg ~ workplace.avg, mobility.active.Vermont)
workplace.active.Vermont.sum <- summary(workplace.active.Vermont.lm)

residential.active.Vermont.lm <- lm(active.avg ~ residential.avg, mobility.active.Vermont)
residential.active.Vermont.sum <- summary(residential.active.Vermont.lm)
#significant results: retail.active.Vermont.sum (***), grocery.active.Vermont.sum (***),parks.active.Vermont.sum (**), transit.active.Vermont.sum (***), workplace.active.Vermont.sum (**), residential.active.Vermont.sum(***)

Iowa

mobility.active.Iowa <- mobility.active %>%
  filter(state == "Iowa")

retail.active.Iowa.lm <- lm(active.avg ~ retail.recreation.avg, mobility.active.Iowa)
retail.active.Iowa.sum <- summary(retail.active.Iowa.lm)

grocery.active.Iowa.lm <- lm(active.avg ~ grocery.pharmacy.avg, mobility.active.Iowa)
grocery.active.Iowa.sum <- summary(grocery.active.Iowa.lm)

parks.active.Iowa.lm <- lm(active.avg ~ parks.avg, mobility.active.Iowa)
parks.active.Iowa.sum <- summary(parks.active.Iowa.lm)

transit.active.Iowa.lm <- lm(active.avg ~ transit.avg, mobility.active.Iowa)
transit.active.Iowa.sum <- summary(transit.active.Iowa.lm)

workplace.active.Iowa.lm <- lm(active.avg ~ workplace.avg, mobility.active.Iowa)
workplace.active.Iowa.sum <- summary(workplace.active.Iowa.lm)

residential.active.Iowa.lm <- lm(active.avg ~ residential.avg, mobility.active.Iowa)
residential.active.Iowa.sum <- summary(residential.active.Iowa.lm)
#significant results: retail.active.Iowa.sum (*), workplace.active.Iowa.sum (***),residential.active.Iowa.sum (*)

Graphs

Cumulative death cases

covid.death.midwest

covid.death.northeast

covid.death.south

covid.death.west

Cumulative confirmed cases

covid.confirm.midwest

covid.confirm.northeast

covid.confirm.south

covid.confirm.west

Testing rates

covid.test.rate.midwest

covid.test.rate.northeast

covid.test.rate.south

covid.test.rate.west

Positive testing rates

covid.incident.rate.midwest 

covid.incident.rate.northeast

covid.incident.rate.south

covid.incident.rate.west

Mean Percent Change of Retail and Recreation

mobility.r.midwest

mobility.r.northeast

mobility.r.south

mobility.r.west

Mean Percent Change of Grocery and Pharmacy

mobility.gp.midwest

mobility.gp.northeast

mobility.gp.south

mobility.gp.west

Mean Percent Change of Parks

mobility.parks.midwest

## Warning: Removed 13 rows containing missing values (geom_point).

mobility.parks.northeast

## Warning: Removed 10 rows containing missing values (geom_point).

mobility.parks.south

## Warning: Removed 12 rows containing missing values (geom_point).

## Warning: Removed 7 row(s) containing missing values (geom_path).

## Warning: Removed 7 row(s) containing missing values (geom_path).

mobility.parks.west

## Warning: Removed 14 rows containing missing values (geom_point).

Mean Percent Change of Transit Stations

mobility.transit.midwest

## Warning: Removed 2 rows containing missing values (geom_point).

mobility.transit.northeast

## Warning: Removed 6 rows containing missing values (geom_point).

mobility.transit.south

## Warning: Removed 6 rows containing missing values (geom_point).

mobility.transit.west

## Warning: Removed 5 rows containing missing values (geom_point).

Mean Percent Change of Workplaces

mobility.work.midwest

mobility.work.northeast

mobility.work.south

mobility.work.west

Mean Percent Change of Residential (%change of people staying at home)

mobility.res.midwest 

## Warning: Removed 1 rows containing missing values (geom_point).

mobility.res.northeast

## Warning: Removed 4 rows containing missing values (geom_point).

mobility.res.south

## Warning: Removed 3 rows containing missing values (geom_point).

mobility.res.west

## Warning: Removed 2 rows containing missing values (geom_point).

2020 unemployment rate

unemploy.mw

unemploy.ne

unemploy.s

unemploy.w

Poverty rate

svi.pov.mw

svi.pov.ne

svi.pov.s

svi.pov.w

## Warning: Removed 1 rows containing missing values (position_stack).

Uninsured rate

svi.uninsured.mw

svi.uninsured.ne

svi.uninsured.s

svi.uninsured.w

Mobile home rate

svi.mobile.mw

svi.mobile.ne

svi.mobile.s

svi.mobile.w

Minority rate

svi.minority.mw

svi.minority.ne

svi.minority.s

svi.minority.w

Disabled rate

svi.disabled.mw

svi.disabled.ne

svi.disabled.s

svi.disabled.w

Elderly rate

svi.old.mw

svi.old.ne

svi.old.s

svi.old.w

Per capita income

svi.pci.mw

svi.pci.ne

svi.pci.s

svi.pci.w

## Warning: Removed 1 rows containing missing values (position_stack).